#-------------------------------------------------------------------------------*
# Natural Resource Rents, Local Taxes, and Government Performance: Evidence from
# Colombia
# By: Luis R. Martinez
# Name: 31-appendix-maps.R
# Description: This do file produces both maps in figure B3.
#-------------------------------------------------------------------------------*
  
cat("\014")

# Define you path to the folder here, making sure to use forward slashes (/):
p <- "XXX"

# Required Packages:
install.packages("sp")
install.packages("maptools")
install.packages("rgeos")
install.packages("rgdal")
install.packages("ggplot2")
install.packages("mapproj")
install.packages("plyr")

library(sp)
library(rgdal)
library(ggplot2)
library(maptools) # used here to read shapefiles
library(rgeos)
library(mapproj)
library(plyr)
library(RColorBrewer)
library(scales)

# Importing official shapefiles:
col <- readShapePoly(file.path(p,"data/raw/13-maps/COL_adm0"))
dpt <- readShapePoly(file.path(p,"data/raw/13-maps/COL_adm1"))
muni <-readShapePoly(file.path(p,"data/raw/13-maps/COL_adm2"))
ven <- readShapePoly(file.path(p,"data/raw/13-maps/VEN_adm0"))
pan <- readShapePoly(file.path(p,"data/raw/13-maps/PAN_adm0"))
ecu <- readShapePoly(file.path(p,"data/raw/13-maps/ECU_adm0"))
per <- readShapePoly(file.path(p,"data/raw/13-maps/PER_adm0"))
bra <- readShapePoly(file.path(p,"data/raw/13-maps/BRA_adm0"))

# Some municipal codes in the department of Vichada need fixing to merge with
# data on cadastral update and oil production.
fixvichada.df <- read.csv(file.path(p,"data/raw/13-maps/fixvichada.csv"))
muni1 <- muni
muni1 <- merge(muni,fixvichada.df, by = "ID_2") 
muni.ids <- muni1[["ID"]]
muni_fix <- gUnionCascaded(muni1, id = muni.ids)

# Define data set fo the maps:
col.df <- fortify(col)
dpt.df <- fortify(dpt)
muni.df <- fortify(muni_fix)
ven.df <- fortify(ven)
pan.df <- fortify(pan)
ecu.df <- fortify(ecu)
per.df <- fortify(per)
bra.df <- fortify(bra)

# Add data con cadastral update and oil-royalties:
data_maps.df <- read.csv(file.path(p,"data/raw/13-maps/data_maps.csv"))
muni_data.df <- merge(muni.df,data_maps.df, by = "id") 
muni_data.df <- join(muni.df,data_maps.df) 

# Generate polygon objects: 
pmuni<-c(geom_polygon(data=muni.df, aes(x=long, y=lat, group=group), colour= "grey", fill="transparent"))
pdpt<-c(geom_polygon(data=dpt.df, aes(x=long, y=lat, group=group), colour= "black", fill="transparent"))
pbog<-c(geom_polygon(data=subset(muni_data.df,id==568), aes(x=long, y=lat, group=group), colour= "black", fill="transparent"))
pmiss<-c(geom_polygon(data=subset(muni_data.df,valid==0), aes(x=long, y=lat, group=group), colour= "grey", fill="#969696"))
pven<-c(geom_polygon(data=ven.df, aes(x=long, y=lat, group=group), colour= "black", fill="#a1d99b"))
ppan<-c(geom_polygon(data=pan.df, aes(x=long, y=lat, group=group), colour= "black", fill="#a1d99b"))
pecu<-c(geom_polygon(data=ecu.df, aes(x=long, y=lat, group=group), colour= "black", fill="#a1d99b"))
pper<-c(geom_polygon(data=per.df, aes(x=long, y=lat, group=group), colour= "black", fill="#a1d99b"))
pbra<-c(geom_polygon(data=bra.df, aes(x=long, y=lat, group=group), colour= "black", fill="#a1d99b"))

# Borders with Colombia:
lpais <- data.frame(
  long = c(-75,-68,-68,-78.1,-77), 
  lat = c(-3,-2,9,9.9,-1),
  group=c(0,0,0,0,0),
  label = c("PERU", "BRAZIL","VENEZUELA","PANAMA","ECUADOR") 
)

#MAP B3-A: Cadastral updates 2006-2010
ggplot(muni_data.df, aes(long,lat, group=group)) + 
  geom_polygon(aes(fill=factor(fecha_up))) + 
  coord_fixed() + 
  pven + ppan + pecu + pper + pbra + pmiss + pmuni + pbog + pdpt+ 
  coord_map(projection="mercator",xlim = c(-79.2,-66.5),ylim = c(-4.5, 13))+
  theme_bw() +guides(fill=guide_legend(title=NULL)) + xlab("longitude") + ylab("latitude")+ 
  scale_fill_manual(values=c("white", "#00ffcc","#0099ff", "#9966ff", "#ff33cc","#ff6600"),breaks=c("0","2006","2007","2008","2009","2010"),labels=c("No update","2006","2007","2008","2009","2010")) 

ggsave(file.path(p,"results/appendix/B-additional/03-figureB3_a.png"), dpi=100)

#MAP B3-B: Oil royalties 2000-2004
ggplot(muni_data.df, aes(long,lat, group=group)) + 
  geom_polygon(aes(fill=factor(roilcut))) + 
  coord_fixed() + 
  pven + ppan + pecu + pper + pbra + pmiss + pmuni + pbog + pdpt+ 
  coord_map(projection="mercator",xlim = c(-79.2,-66.5),ylim = c(-4.5, 13))+
  theme_bw() +guides(fill=guide_legend(title=NULL)) + xlab("longitude") + ylab("latitude")+ 
  scale_fill_manual(values=c("white", "#ffff66","#ffcc66", "#ff9933", "#ff3300","#990000"),breaks=c("0","1","2","3","4","5"),labels=c("No oil","Quintile 1","Quintile 2","Quintile 3","Quintile 4","Quintile 5")) 

ggsave(file.path(p,"results/appendix/B-additional/03-figureB3_b.png"), dpi=100)

